## DGP for simulation studies
dgp <- function (coef.ps, N, formula_y, formula_ps) {
	covariate <- rnorm(n = N * 4, mean = rep(0, 4), sd = rep(1, 4))
	covariate <- matrix(covariate, nrow = N, ncol = 4)
	colnames(covariate) <- paste0("x", 1:4)
	data <- data.frame(covariate) %>%
					mutate(x001 = exp(x1 / 2),
							   x002 = x2 / (1 + exp(x1)) + 10,
							   x003 = (x1 * x3 / 25 + 0.6)^3,
							   x004 = (x2 + x4 + 20)^2)
	assign.ps <- dgp_ps(formula = formula_ps, 
								  		coef = coef.ps,
								  		data = data)
	data$treat <- assign.ps$treat
	data$tps <- assign.ps$ps
	data <- dgp0(formula = formula_y,
					  	 data = data) %>%
					mutate(x1 = scale(x1),
								 x2 = scale(x2),
								 x3 = scale(x3),
								 x4 = scale(x4),
								 x001 = scale(x001),
								 x002 = scale(x002),
								 x003 = scale(x003),
								 x004 = scale(x004))
	data
}

## Data generating functions (treatment assignment)
dgp_ps <- function (formula, coef, data) {
	N <- nrow(data)
  model_dgp <- stats::model.frame(formula, data = data)
  X <- as.matrix(stats::model.matrix(formula, model_dgp))
	y_star <- logistic(X %*% coef)
	treat <- rbinom(n = rep(1, N), size = rep(1, N), prob = y_star)
	list(treat = treat, ps = y_star)
}

## Data generating functions (outcome)
dgp0 <- function (formula, data) {
	coef.y_lin1 <- c(210, 0, 27.4, 13.7, 13.7, 13.7) # for linear1
	coef.y_lin2 <- c(210, 0, 13.7, 27.4, 27.4, 27.4) # for linear2
	coef.y_quadexp <- c(210, 0, 27.4, 27.4, 27.4, 27.4) # for quad1, quad2, exp1, and exp2
	N <- nrow(data)
  model_dgp <- stats::model.frame(formula, data = data)
  X <- as.matrix(stats::model.matrix(formula, model_dgp))

	data$ylin1 <- c(X %*% coef.y_lin1) + rnorm(N, mean = 0, sd = 1)
	data$ylin2 <- c(X %*% coef.y_lin2) + rnorm(N, mean = 0, sd = 1)

  Xquad1 <- X
  Xquad1[Xquad1 < 0] <- 0
	data$yquad1 <- c(Xquad1^2 %*% coef.y_quadexp) + rnorm(N, mean = 0, sd = 1)
  
  Xquad2 <- X
  Xquad2[Xquad2 > 0 & Xquad2 != 1] <- 0
	data$yquad2 <- c(Xquad2^2 %*% coef.y_quadexp) + rnorm(N, mean = 0, sd = 1)
  
  Xexp1 <- X
	Xexp1[, -1] <- exp(Xexp1[, -1])
	data$yexp1 <- c(Xexp1 %*% coef.y_quadexp) + rnorm(N, mean = 0, sd = 1)
  
  Xexp2 <- X
	Xexp2[, -1] <- exp(-Xexp2[, -1])
	data$yexp2 <- c(Xexp2 %*% coef.y_quadexp) + rnorm(N, mean = 0, sd = 1)

	data
}

## Logistic function
logistic <- function (x) {exp(x) / (1 + exp(x))}

## Simulation main functions
sim_main <- function (B, save, seed, lambda200, lambda1000) {
  res200_res <- sim_main_internal(N = 200, B = B, lambda = lambda200, seed = seed)
  print("start 1000")
  res1000_res <- sim_main_internal(N = 1000, B = B, lambda = lambda1000, seed = seed)
  res200 <- res200_res$result
  res1000 <- res1000_res$result
  if (save == TRUE) {
  	sim_tbl1(result200 = res200$lin1, result1000 = res1000$lin1, type = "linear1",
  					 caption = "Linear outcome model 1")
  	sim_tbl1(result200 = res200$lin2, result1000 = res1000$lin2, type = "linear2",
  					 caption = "Linear outcome model 2")
  	sim_tbl1(result200 = res200$quad1, result1000 = res1000$quad1, type = "quad1",
  					 caption = "Quadratic outcome model 1")
  	sim_tbl1(result200 = res200$quad2, result1000 = res1000$quad2, type = "quad2",
  					 caption = "Quadratic outcome model 2")
  	sim_tbl1(result200 = res200$exp1, result1000 = res1000$exp1, type = "exp1",
  					 caption = "Exponential outcome model 1")
  	sim_tbl1(result200 = res200$exp2, result1000 = res1000$exp2, type = "exp2",
  					 caption = "Exponential outcome model 2")
  }
  result <- list(res200 = res200, res1000 = res1000)
  notconverged <- list(notconverged200 = res200_res$notconverged,
  										 notconverged1000 = res1000_res$notconverged)
  list(result = result, notconverged = notconverged)
}

## Simulation result table 1
sim_tbl1 <- function (result200, result1000, type, caption) {
  res200s <- result200 %>%
  					 select(-c(N)) %>%
  					 mutate(Bias = round(Bias, digits = 2), 
  					 				RMSE = round(RMSE, digits = 2)) %>%
  					 rename(Bias200 = Bias, RMSE200 = RMSE)
  res1000s <- result1000 %>%
  					 	select(Bias, RMSE) %>%
  						mutate(Bias = round(Bias, digits = 2), 
  					 				 RMSE = round(RMSE, digits = 2)) %>%
  						rename(Bias1000 = Bias, RMSE1000 = RMSE)

 	res1 <- data.frame(res200s, res1000s) %>%
 					filter(psmodel == "correct" & posneg == "negative") %>%
 					select(Estimator, Bias200, RMSE200, Bias1000, RMSE1000)
 	res1 <- data.frame(res1[, 2:3], NA, res1[, 4:5], NA)

 	res2 <- data.frame(res200s, res1000s) %>%
 					filter(psmodel != "correct" & posneg == "negative") %>%
 					select(Estimator, Bias200, RMSE200, Bias1000, RMSE1000)
 	res2 <- data.frame(res2[, 2:3], NA, res2[, 4:5], NA)

 	res3 <- data.frame(res200s, res1000s) %>%
 					filter(psmodel == "correct" & posneg == "positive") %>%
 					select(Estimator, Bias200, RMSE200, Bias1000, RMSE1000)
 	res3 <- data.frame(res3[, 2:3], NA, res3[, 4:5])

 	res4 <- data.frame(res200s, res1000s) %>%
 					filter(psmodel != "correct" & posneg == "positive") %>%
 					select(Estimator, Bias200, RMSE200, Bias1000, RMSE1000)
 	res4 <- data.frame(res4[, 2:3], NA, res4[, 4:5])

 	result1 <- rbind(res1, res2)
 	result2 <- rbind(res3, res4)
 	result <- data.frame(result1, NA, result2)

	clabels <- c(rep(c("Bias", "RMSE", ""), 2), rep(c("", "Bias", "RMSE"), 2))
	caption <- paste("Simulation results:", caption)
	rownames <- res200s$Estimator[c(1:14, 29:42)]
	rownames[rownames == "nDBW"] <- "\\textbf{nDBW}"
	rownames[rownames == "nDBW DR"] <- "\\textbf{nDBW DR}"
	note <- 
		"\\parbox{0.99\\textwidth}
		{Notes: This simulation compares the performance of various methods 
		for estimating propensity scores and (inverse probability) weights 
		by investigating combinations of six versions of the true outcome model 
		(linear~1, linear~2, quadratic~1, quadratic~2, exponential~1, and exponential~2)
		and two versions of coefficients for the true propensity score model (type~A and B)
		with the two different numbers of observations ($n = 200$ and $n = 1000$).
		For each estimation method, I use two propensity score model specifications 
		(correct and misspecified) and report the bias and RMSE for each in the table.}"

	filename <- "tab_i1"
	if (type == "linear2") {
		filename <- "tab_i2"
	} else if (type == "quad1") {
		filename <- "tab_i3"
	} else if (type == "quad2") {
		filename <- "tab_i4"
	} else if (type == "exp1") {
		filename <- "tab_i5"
	} else if (type == "exp2") {
		filename <- "tab_i6"
	}

	## LaTeX table
	latex(object = result,
	      title = "",
	      file = paste0("tables/", filename, ".tex"),
	      #file = paste0("tables/", type, ".tex"),
	      label = paste0("tb_", type),
	      caption = caption,
	      insert.bottom = note,
	      first.hline.double = FALSE,
	      rowname = rownames,
	      cgroup = c("$n = 200$", "", "$n = 1000$", "",
	      					 "$n = 200$", "", "$n = 1000$"),
	      n.cgroup = c(2, 1, 2, 2, 2, 1, 2),
	      cgroupTexCmd = "",
	      colheads = clabels,
	      rgroup = c("Correct PS model", "Misspecified PS model"),
	      n.rgroup = c(rep(nrow(res1), 2)),
	      longtable = FALSE,
	      center = "centering")
}

## Simulation result table 2
sim_tbl2 <- function (result1, result2, type) {
	tbl2 <- c("nDBW DR", "MLE DR", "CBPS DR/BRDR", "CBPS DR", "Calibrated weighting DR", 
						"Entropy balancing DR", "True propensity score DR~~", "Unweighted")
  res200s1 <- result1 %>%
  					  filter(N == 200) %>%
  					  select(-c(N)) %>%
  					  mutate(Bias = round(Bias, digits = 2), 
  					 				 RMSE = round(RMSE, digits = 2)) %>%
  					  rename(Bias200 = Bias, RMSE200 = RMSE) %>%
  					  filter(Estimator %in% tbl2)
  res1000s1 <- result1 %>%
  						 filter(N == 1000) %>%
  					   filter(Estimator %in% tbl2) %>%
  						 select(Bias, RMSE) %>%
  						 mutate(Bias = round(Bias, digits = 2), 
  						 			  RMSE = round(RMSE, digits = 2)) %>%
  						 rename(Bias1000 = Bias, RMSE1000 = RMSE)
  res200s2 <- result2 %>%
  					  filter(N == 200) %>%
  					  select(-c(N)) %>%
  					  mutate(Bias = round(Bias, digits = 2), 
  					 				 RMSE = round(RMSE, digits = 2)) %>%
  					  rename(Bias200 = Bias, RMSE200 = RMSE) %>%
  					  filter(Estimator %in% tbl2)
  res1000s2 <- result2 %>%
  						 filter(N == 1000) %>%
  					   filter(Estimator %in% tbl2) %>%
  						 select(Bias, RMSE) %>%
  						 mutate(Bias = round(Bias, digits = 2), 
  						 			  RMSE = round(RMSE, digits = 2)) %>%
  						 rename(Bias1000 = Bias, RMSE1000 = RMSE)

  res11 <- data.frame(res200s1, res1000s1) %>%
  				 filter(psmodel == "correct" & posneg == "negative") %>%
  				 select(Estimator, Bias200, RMSE200, Bias1000, RMSE1000)
  res11 <- data.frame(res11[, 2:3], NA, res11[, 4:5], NA)

  res12 <- data.frame(res200s1, res1000s1) %>%
  				 filter(psmodel != "correct" & posneg == "negative") %>%
  				 select(Estimator, Bias200, RMSE200, Bias1000, RMSE1000)
  res12 <- data.frame(res12[, 2:3], NA, res12[, 4:5], NA)

  res13 <- data.frame(res200s1, res1000s1) %>%
  				 filter(psmodel == "correct" & posneg == "positive") %>%
  				 select(Estimator, Bias200, RMSE200, Bias1000, RMSE1000)
  res13 <- data.frame(res13[, 2:3], NA, res13[, 4:5])

  res14 <- data.frame(res200s1, res1000s1) %>%
  				 filter(psmodel != "correct" & posneg == "positive") %>%
  				 select(Estimator, Bias200, RMSE200, Bias1000, RMSE1000)
  res14 <- data.frame(res14[, 2:3], NA, res14[, 4:5])

  result_11 <- rbind(res11, res12)
  result_12 <- rbind(res13, res14)
  result_1 <- data.frame(result_11, NA, result_12)

  res21 <- data.frame(res200s2, res1000s2) %>%
  				 filter(psmodel == "correct" & posneg == "negative") %>%
  				 select(Estimator, Bias200, RMSE200, Bias1000, RMSE1000)
  res21 <- data.frame(res21[, 2:3], NA, res21[, 4:5], NA)

  res22 <- data.frame(res200s2, res1000s2) %>%
  				 filter(psmodel != "correct" & posneg == "negative") %>%
  				 select(Estimator, Bias200, RMSE200, Bias1000, RMSE1000)
  res22 <- data.frame(res22[, 2:3], NA, res22[, 4:5], NA)

  res23 <- data.frame(res200s2, res1000s2) %>%
  				 filter(psmodel == "correct" & posneg == "positive") %>%
  				 select(Estimator, Bias200, RMSE200, Bias1000, RMSE1000)
  res23 <- data.frame(res23[, 2:3], NA, res23[, 4:5])

  res24 <- data.frame(res200s2, res1000s2) %>%
  				 filter(psmodel != "correct" & posneg == "positive") %>%
  				 select(Estimator, Bias200, RMSE200, Bias1000, RMSE1000)
  res24 <- data.frame(res24[, 2:3], NA, res24[, 4:5])

  result_21 <- rbind(res21, res22)
  result_22 <- rbind(res23, res24)
  result_2 <- data.frame(result_21, NA, result_22)

  result <- rbind(result_1, result_2)

	clabels <- c(rep(c("Bias", "RMSE", ""), 2), rep(c("", "Bias", "RMSE"), 2))
	caption <- paste("Simulation results:", type, "outcome models")
	rownames <- res200s1$Estimator[c(1:7, 15:21, 8:14, 22:28)]
	rownames[rownames == "nDBW DR"] <- "\\textbf{nDBW DR}"
	note <- 
	paste0("\\parbox{0.99\\textwidth}
	{Notes: This simulation compares the performance of various methods 
	for estimating propensity scores and (inverse probability) weights 
	by investigating combinations of two versions of the true outcome model 
	(", type, "~1 and 2)
	and two versions of coefficients for the propensity score model (type~A and B)
	with the two different numbers of observations ($n = 200$ and $n = 1000$).
	For each estimation method, I use two propensity score model specifications 
	(correct and misspecified) and report the bias and RMSE for each in the table.}")
	rgroups <- c(paste0("\\multicolumn{14}{l}{\\textbf{", type, " outcome model 1: correct PS model}"),
							 paste0("\\multicolumn{14}{l}{\\textbf{", type, " outcome model 1: misspecified PS model}"),
							 paste0("\\multicolumn{14}{l}{\\textbf{", type, " outcome model 2: correct PS model}"),
							 paste0("\\multicolumn{14}{l}{\\textbf{", type, " outcome model 2: misspecified PS model}"))

	filename <- "tab_2"
	if (type == "Quadratic") {
		filename <- "tab_3"
	} else if (type == "Exponential") {
		filename <- "tab_4"
	}

	## LaTeX table 2
	latex(object = result,
		    title = "",
		    file = paste0("tables/", filename, ".tex"),
		    #file = paste0("tables/", type, "_2.tex"),
		    label = paste0("tb_", type, "_2"),
		    caption = caption,
		    insert.bottom = note,
		    first.hline.double = FALSE,
		    rowname = rownames,
		    cgroup = c("$n = 200$", "", "$n = 1000$", "",
		      				 "$n = 200$", "", "$n = 1000$"),
		    n.cgroup = c(2, 1, 2, 2, 2, 1, 2),
		    cgroupTexCmd = "",
		    colheads = clabels,
		    rgroup = rgroups,
		    n.rgroup = c(rep(7, 4)),
		    rgroupTexCmd = "",
		    longtable = FALSE,
		    center = "centering")
}


sim_main_internal <- function (N, B, lambda, seed) {
	result <- future_map2(.x = rep(c("correct", "missp"), each = 2), 
												.y = rep(c("positive", "negative"), times = 2), 
												.f = sim_dbw, 
												N = N, B = B, lambda = lambda,
												.options = furrr_options(seed = seed))
	corpos <- result[[1]]
	corneg <- result[[2]]
	mispos <- result[[3]]
	misneg <- result[[4]]
  lin1 <- data.frame(N = N, corpos$summary$lin1) %>%
					bind_rows(data.frame(N = N, corneg$summary$lin1)) %>%
					bind_rows(data.frame(N = N, mispos$summary$lin1)) %>%
					bind_rows(data.frame(N = N, misneg$summary$lin1))
  lin2 <- data.frame(N = N, corpos$summary$lin2) %>%
					bind_rows(data.frame(N = N, corneg$summary$lin2)) %>%
					bind_rows(data.frame(N = N, mispos$summary$lin2)) %>%
					bind_rows(data.frame(N = N, misneg$summary$lin2))
  quad1 <- data.frame(N = N, corpos$summary$quad1) %>%
					 bind_rows(data.frame(N = N, corneg$summary$quad1)) %>%
					 bind_rows(data.frame(N = N, mispos$summary$quad1)) %>%
					 bind_rows(data.frame(N = N, misneg$summary$quad1))
  quad2 <- data.frame(N = N, corpos$summary$quad2) %>%
					 bind_rows(data.frame(N = N, corneg$summary$quad2)) %>%
					 bind_rows(data.frame(N = N, mispos$summary$quad2)) %>%
					 bind_rows(data.frame(N = N, misneg$summary$quad2))
  exp1 <- data.frame(N = N, corpos$summary$exp1) %>%
					bind_rows(data.frame(N = N, corneg$summary$exp1)) %>%
					bind_rows(data.frame(N = N, mispos$summary$exp1)) %>%
					bind_rows(data.frame(N = N, misneg$summary$exp1))
  exp2 <- data.frame(N = N, corpos$summary$exp2) %>%
					bind_rows(data.frame(N = N, corneg$summary$exp2)) %>%
					bind_rows(data.frame(N = N, mispos$summary$exp2)) %>%
					bind_rows(data.frame(N = N, misneg$summary$exp2))
	notconverged <- list(corpos = corpos$notconverged,
											 corneg = corneg$notconverged,
											 mispos = mispos$notconverged,
											 misneg = misneg$notconverged)
	result <- list(lin1 = lin1, lin2 = lin2, quad1 = quad1, quad2 = quad2, exp1 = exp1, exp2 = exp2)
	list(result = result, notconverged = notconverged)
}


## Simulation internal functions
sim_dbw <- function (psmodel, coef.ps_pn, N, B, lambda) {
	truelin <- 210 # for linear1 and linear2
	truequad <- 210 + 27.4 * 4 / 2 # for quad1 and quad2
	trueexp <- 210 + 27.4 * 4 * exp(1 / 2) # for exp1 and exp2
	formula_trueps <- as.formula("~ x1 + x2 + x3 + x4")
	formula_truey <- as.formula("~ treat + x1 + x2 + x3 + x4")
	formula_y <- as.formula("ylin1 ~ treat + x001 + x002 + x003 + x004")
	if (coef.ps_pn == "positive") {
		coef.ps <- c(0.0, 1.0, -0.5, 0.25, 0.1)
	} else {
		coef.ps <- -c(0.0, 1.0, -0.5, 0.25, 0.1)
	}
	if (psmodel == "correct") {
		formula_ps <- as.formula("treat ~ x1 + x2 + x3 + x4")
	} else {
		formula_ps <- as.formula("treat ~ x001 + x002 + x003 + x004")
	}
	formula_y_temp <- as.formula("ylin1 ~ x001 + x002 + x003 + x004")
  res_ndbw <- list()
  res_mle <- list()
  res_cbps <- list()
  res_cb <- list()
  res_ebal <- list()
  res_tps <- list()
  res_unw <- list()
  converged <- numeric(B)
	for (b in 1:B) {
    df <- dgp(coef.ps = coef.ps, N = N, formula_y = formula_truey, formula_ps = formula_trueps)
    result_ndbw <- dbw(formula_y = formula_y_temp, 
                   		 formula_ps = formula_ps, 
                  		 estimand = "AO",
                  		 method = "dbw",
                  		 method_y = "wls",
                  		 data = df,
                  		 normalize = TRUE,
                  		 vcov = FALSE,
                  		 lambda = lambda,
                  		 tol = 1e-7,
                  		 init_lambda = 0.05)
    converged[b] <- result_ndbw$converged
    result_mle <- dbw(formula_y = formula_y_temp, 
                   		formula_ps = formula_ps, 
                  		estimand = "AO",
                  		method = "mle",
                  		method_y = "wls",
                  		data = df,
                  		normalize = FALSE,
                  		vcov = FALSE,
                  		lambda = 0)
    result_cbps <- dbw(formula_y = formula_y_temp, 
                   		 formula_ps = formula_ps, 
                  		 estimand = "ATEcombined",
                  		 method = "cb",
                  		 method_y = "wls",
                  		 data = df,
                  		 normalize = TRUE,
                  		 vcov = FALSE,
                  		 lambda = 0)
    result_cb <- dbw(formula_y = formula_y_temp, 
                   	 formula_ps = formula_ps, 
                  	 estimand = "AO",
                 		 method = "cb",
                 		 method_y = "wls",
                 		 data = df,
                 		 normalize = TRUE,
                  	 vcov = FALSE,
                 		 lambda = 0)
    result_ebal <- dbw(formula_y = formula_y_temp, 
                 	  	 formula_ps = formula_ps, 
                  		 estimand = "AO",
                 			 method = "eb",
                	 		 method_y = "wls",
                 			 data = df,
                 			 normalize = TRUE,
                  		 vcov = FALSE,
                 			 lambda = 0)
    response <- result_ndbw$response
  	weights <- rep(1, N)
		res_ndbw[[b]] <- sim_dbw_y(estimator = c("nDBW", "nDBWDR"), df = df, 
															 response = response, weights = result_ndbw$est_weights, N = N, 
															 est_weights = result_ndbw$est_weights)
		res_mle[[b]] <- sim_dbw_y(estimator = c("MLE", "MLEDR"), df = df, 
															response = response, weights = weights, N = N, 
															est_weights = result_mle$est_weights)
		res_cbps[[b]] <- sim_dbw_y(estimator = c("CBPS", "CBPSDR"), df = df, 
															 response = response, weights = weights, N = N, 
															 est_weights = result_cbps$est_weights)
		res_cb[[b]] <- sim_dbw_y(estimator = c("CBW", "CBWDR"), df = df, 
														 response = response, weights = weights, N = N, 
														 est_weights = result_cb$est_weights)
		res_ebal[[b]] <- sim_dbw_y(estimator = c("EB", "EBDR"), df = df, 
															 response = response, weights = weights, N = N, 
															 est_weights = result_ebal$est_weights)
		res_tps[[b]] <- sim_dbw_y(estimator = c("TPS", "TPSDR"), df = df, 
															response = response, weights = weights, N = N, 
															est_weights = 1 / df$tps / N)
		res_unw[[b]] <- sim_dbw_y(estimator = c("Unw", "UnwDR"), df = df, 
															response = response, weights = weights, N = N, 
															est_weights = rep(1 / sum(response), N))
	}
	res_ndbw <- bind_rows(res_ndbw)
	res_mle <- bind_rows(res_mle)
	res_cbps <- bind_rows(res_cbps)
	res_cb <- bind_rows(res_cb)
	res_ebal <- bind_rows(res_ebal)
	res_tps <- bind_rows(res_tps)
	res_unw <- bind_rows(res_unw)
	result <- rbind(res_ndbw, res_mle, res_cbps, res_cb, res_ebal, res_tps, res_unw)
	lin1 <- subset(result, type == "lin1")
	lin2 <- subset(result, type == "lin2")
	quad1 <- subset(result, type == "quad1")
	quad2 <- subset(result, type == "quad2")
	exp1 <- subset(result, type == "exp1")
	exp2 <- subset(result, type == "exp2")
	lin1 <- est_summarise(result = lin1, type = "linear1", true = truelin, 
												coef.ps_pn = coef.ps_pn, psmodel = psmodel)
	lin2 <- est_summarise(result = lin2, type = "linear2", true = truelin, 
												coef.ps_pn = coef.ps_pn, psmodel = psmodel)
	quad1 <- est_summarise(result = quad1, type = "quad1", true = truequad, 
												 coef.ps_pn = coef.ps_pn, psmodel = psmodel)
	quad2 <- est_summarise(result = quad2, type = "quad2", true = truequad, 
												 coef.ps_pn = coef.ps_pn, psmodel = psmodel)
	exp1 <- est_summarise(result = exp1, type = "exp1", true = trueexp, 
												coef.ps_pn = coef.ps_pn, psmodel = psmodel)
	exp2 <- est_summarise(result = exp2, type = "exp2", true = trueexp, 
												coef.ps_pn = coef.ps_pn, psmodel = psmodel)
	est_summary <- list(lin1 = lin1, lin2 = lin2, quad1 = quad1, quad2 = quad2, exp1 = exp1, exp2 = exp2)
	list(summary = est_summary, notconverged = which(converged == 0))
}

## Simulation function (internal)
sim_dbw_y <- function (estimator, df, response, weights, N, est_weights) {
	formula_y_lin1 <- stats::as.formula(ylin1 ~ x001 + x002 + x003 + x004)
	formula_y_lin2 <- stats::as.formula(ylin2 ~ x001 + x002 + x003 + x004)
	formula_y_quad1 <- stats::as.formula(yquad1 ~ x001 + x002 + x003 + x004)
	formula_y_quad2 <- stats::as.formula(yquad2 ~ x001 + x002 + x003 + x004)
	formula_y_exp1 <- stats::as.formula(yexp1 ~ x001 + x002 + x003 + x004)
	formula_y_exp2 <- stats::as.formula(yexp2 ~ x001 + x002 + x003 + x004)
  result_y_lin1 <- dbw_y_wls(formula_y = formula_y_lin1, 
  													 df = df[response == 1, ], 
  													 fulldf = df, 
  													 est_weights = weights[response == 1],
  													 response = response)
  result_y_lin2 <- dbw_y_wls(formula_y = formula_y_lin2, 
  													 df = df[response == 1, ], 
  													 fulldf = df, 
  													 est_weights = weights[response == 1],
  													 response = response)
  result_y_quad1 <- dbw_y_wls(formula_y = formula_y_quad1, 
  													 df = df[response == 1, ], 
  													 fulldf = df, 
  													 est_weights = weights[response == 1],
  													 response = response)
  result_y_quad2 <- dbw_y_wls(formula_y = formula_y_quad2, 
  													 df = df[response == 1, ], 
  													 fulldf = df, 
  													 est_weights = weights[response == 1],
  													 response = response)
  result_y_exp1 <- dbw_y_wls(formula_y = formula_y_exp1, 
  													 df = df[response == 1, ], 
  													 fulldf = df, 
  													 est_weights = weights[response == 1],
  													 response = response)
  result_y_exp2 <- dbw_y_wls(formula_y = formula_y_exp2, 
  													 df = df[response == 1, ], 
  													 fulldf = df, 
  													 est_weights = weights[response == 1],
  													 response = response)
  est_lin1 <- sum(response * result_y_lin1$y * est_weights)
  est_dr_lin1 <- sum(response * (result_y_lin1$y - result_y_lin1$pred) * est_weights + 
  										result_y_lin1$pred / N)
  est_lin2 <- sum(response * result_y_lin2$y * est_weights)
  est_dr_lin2 <- sum(response * (result_y_lin2$y - result_y_lin2$pred) * est_weights + 
  										result_y_lin2$pred / N)
  est_quad1 <- sum(response * result_y_quad1$y * est_weights)
  est_dr_quad1 <- sum(response * (result_y_quad1$y - result_y_quad1$pred) * est_weights + 
  										result_y_quad1$pred / N)
  est_quad2 <- sum(response * result_y_quad2$y * est_weights)
  est_dr_quad2 <- sum(response * (result_y_quad2$y - result_y_quad2$pred) * est_weights + 
  										result_y_quad2$pred / N)
  est_exp1 <- sum(response * result_y_exp1$y * est_weights)
  est_dr_exp1 <- sum(response * (result_y_exp1$y - result_y_exp1$pred) * est_weights + 
  										result_y_exp1$pred / N)
  est_exp2 <- sum(response * result_y_exp2$y * est_weights)
  est_dr_exp2 <- sum(response * (result_y_exp2$y - result_y_exp2$pred) * est_weights + 
  										result_y_exp2$pred / N)
  data.frame(estimator = rep(estimator, 6),
  					 type = rep(c("lin1", "lin2", "quad1", "quad2", "exp1", "exp2"), each = 2),
  					 est = c(est_lin1, est_dr_lin1, est_lin2, est_dr_lin2, 
  					 				 est_quad1, est_dr_quad1, est_quad2, est_dr_quad2, 
  					 				 est_exp1, est_dr_exp1, est_exp2, est_dr_exp2))
}

est_summarise <- function (result, type, true, coef.ps_pn, psmodel) {
	est_summary <- df_rmse(estimator = "nDBW", true = true, 
												 est = subset(result, estimator == "nDBW")$est) %>%
								 bind_rows(df_rmse(estimator = "MLE", true = true, 
								 									 est = subset(result, estimator == "MLE")$est)) %>%
								 bind_rows(df_rmse(estimator = "CBPS", true = true, 
								 									 est = subset(result, estimator == "CBPS")$est)) %>%
								 bind_rows(df_rmse(estimator = "Calibrated weighting", true = true, 
								 									 est = subset(result, estimator == "CBW")$est)) %>%
								 bind_rows(df_rmse(estimator = "Entropy balancing", true = true, 
								 									 est = subset(result, estimator == "EB")$est)) %>%
								 bind_rows(df_rmse(estimator = "True propensity score", true = true, 
								 									 est = subset(result, estimator == "TPS")$est)) %>%
								 bind_rows(df_rmse(estimator = "Unweighted", true = true,
								 									 est = subset(result, estimator == "Unw")$est)) %>%
								 bind_rows(df_rmse(estimator = "nDBW DR", true = true,
								 									 est = subset(result, estimator == "nDBWDR")$est)) %>%
								 bind_rows(df_rmse(estimator = "MLE DR", true = true,
								 									 est = subset(result, estimator == "MLEDR")$est)) %>%
								 bind_rows(df_rmse(estimator = "CBPS DR/BRDR", true = true,
								 									 est = subset(result, estimator == "CBPSDR")$est)) %>%
								 bind_rows(df_rmse(estimator = "Calibrated weighting DR", true = true,
								 									 est = subset(result, estimator == "CBWDR")$est)) %>%
								 bind_rows(df_rmse(estimator = "Entropy balancing DR", true = true,
								 									 est = subset(result, estimator == "EBDR")$est)) %>%
								 bind_rows(df_rmse(estimator = "True propensity score DR~~", true = true,
								 									 est = subset(result, estimator == "TPSDR")$est)) %>%
								 bind_rows(df_rmse(estimator = "Imputation", true = true,
								 									 est = subset(result, estimator == "UnwDR")$est))
	est_summary$type <- type
	est_summary$posneg <- coef.ps_pn
	est_summary$psmodel <- psmodel
  est_summary$Estimator[est_summary$Estimator == "CBPS DR/BRDR" & est_summary$psmodel == "correct"] <- "CBPS DR"
	est_summary
}

df_rmse <- function (estimator, true, est) {
	data.frame(Estimator = estimator, Bias = mean(est - true), RMSE = sqrt(mean((est - true)^2)))
}

dbw_y_wls <- function (formula_y, df, fulldf, est_weights, response) {
  df <- data.frame(df, est_weights = est_weights)
  result <- stats::lm(formula = formula_y, data = df, 
                      weights = est_weights, y = TRUE)
  pred <- stats::predict(object = result, newdata = fulldf, type = "response")
  y <- rep(0, length(response))
  y[response == 1] <- result$y
  list(pred = pred, y = y)
}
